Monocle2 踩坑教程（1）

Original 周运来单细胞天地 2022-06-07

分享是一种态度

拟时分析

拟时（pseudotime）分析，又称细胞轨迹（cell trajectory）分析，通过拟时分析可以推断出发育过程细胞的分化轨迹或细胞亚型的演化过程，在发育相关研究中使用频率较高。主要基于关键基因的表达模式，在拟时间中对单个细胞进行排序，模拟出时间发育过程的动态变化。

既然讲到所谓的拟时分析不过是一种排序而已，那么我们就要知道排序的基本要素：

对什么排序
如何判断先后顺序
如何寻找分支点（如果分支的话）

早在20世纪30年代，数量生态学领域前苏联学者Ranensky就提出了排序的概念，并发展了一个简单的排序方法（见Sobolev和Utekhin 1973），但只限于在前苏联传播（Greig-Smith 1980），Ramensky当时应用一个或两个环境因子梯度去排列植物群落，他用的名词是德文“ordnung”。经过多年的发展，排序技术依然是一种在低维空间排布高维数据的降维技术。当我们讲到排序的时候，离不开降维；讲到降维的时候，离不开特征提取（或者选择）。

那么对应到拟时分析的描述中：1)关键基因就是特征选择的结果，2)拟时间就是排序空间，3)排序就是细胞的演化轨迹。所有的拟时分析都离不开这三点。具体地：

Monocle引入了在伪时间(拟时间)内对单个细胞排序的策略，利用单个细胞的非同步进程，将它们置于与细胞分化等生物学过程相对应的轨迹上。Monocle利用先进的机器学习技术(反向图嵌入)从单细胞数据中学习显式的主图来对细胞进行排序，该方法能够可靠、准确地解决复杂的生物学过程。

Monocle2 可以做什么？

Clustering, classifying, and counting cells
Constructing single-cell trajectories.
Differential expression analysis

当然我们关心的是第二个功能了，但是不防也看一下它的其他功能。算法原理可以参看文章Reversed graph embedding resolves complex single-cell developmental trajectories(https://links.jianshu.com/go?to=https%3A%2F%2Fwww.biorxiv.org%2Fcontent%2Fearly%2F2017%2F02%2F21%2F110668.full.pdf)

Monocle 2通过反向图嵌入学习单细胞轨迹的方法。每个细胞都可以表示为高维空间中的一个点，其中每个维对应于有序基因的表达水平。高维数据首先通过PCA(默认方法)、扩散图(diffusion maps)等几种降维方法中的任何一种投影到较低的维空间。然后Monocle 2在自动选择的数据中心集上构造一个生成树。然后，该算法细胞移动到树中最近的顶点，不断更新顶点的位置以“适合”细胞，学习新的生成树，并迭代地继续这个过程，直到树和细胞的位置收敛为止。在整个过程中，Monocle 2保持了高维空间和低维空间之间的可逆映射，从而既学习了轨迹，又降低了数据的维数。一旦Monocle 2习得这棵树，并选择一个tip作为“根”。每个细胞的伪时间计算为其沿树到根的测地距离（Geodesic Distance，在图论中，Geodesic Distance 就是图中两节点的最短路径的距离(https://links.jianshu.com/go?to=http%3A%2F%2Flemonc.me%2Faverage-geodesic-distance.html)），并根据主图自动分配其分支。

主要过程如下，最难理解的应该是反转图嵌入，对算法感兴趣的同学可以研究一下流行学习的知识。

dpFeature selection， identifying genes that define biological process 确定和生物过程有关的基因 selecting the genes differentially expressed between clusters of cells identified with tSNE dimension reduction followed by density peak clustering选择差异表达基因；
反向图嵌入（GRE）：默认DDRTree。它是一种流形学习算法，旨在将主要图形嵌入高维单细胞RNA-seq数据中，找到高维基因表达空间与低维空间映射，同时在这个缩小的空间中学习图的结构（K个质心）来解决这个问题；
伪时间分配和分支计算：利用DDRTree构建MST，将细胞投影到MST上，MST递归计算伪时间；

算法层面适可而止。

seurat2monocle

 1library(monocle)
 2#<-importCDS(pbmc,import_all = TRUE)
 3#Load Seurat object
 4pbmc <- readRDS('D:\\Users\\Administrator\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19pbmc_tutorial.rds')
 5
 6#Extract data, phenotype data, and feature data from the SeuratObject
 7data <- as(as.matrix(pbmc@assays$RNA@counts), 'sparseMatrix')
 8pd <- new('AnnotatedDataFrame', data = pbmc@meta.data)
 9fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
10fd <- new('AnnotatedDataFrame', data = fData)
11
12#Construct monocle cds
13monocle_cds <- newCellDataSet(data,
14                              phenoData = pd,
15                              featureData = fd,
16                              lowerDetectionLimit = 0.5,
17                              expressionFamily = negbinomial.size())

如果你只有count矩阵可以这么读入：

1#do not run
2HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
3HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
4HSMM_gene_annotation <- read.delim("gene_annotations.txt")
5
6pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
7fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
8HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
9    phenoData = pd, featureData = fd)

这里重点是理解CellDataSet 的数据结构。

exprs, a numeric matrix of expression values, where rows are genes, and columns are cells
phenoData, an AnnotatedDataFrame object, where rows are cells, and columns are cell attributes** (such as cell type, culture condition, day captured, etc.)
featureData, an AnnotatedDataFrame object, where rows are features (e.g. genes), and columns are gene attributes, such as biotype, gc content, etc.

官网给出许多其他建议以及读入数据的注意事项，需要注意的是，如果是UMI数据，不应该在创建CellDataSet之前自己对其进行标准化。也不应该试图将UMI计数转换为相对丰度(通过将其转换为FPKM/TPM数据)。Monocle将在内部执行所有需要的标准化步骤。自己将其正常化可能会破坏Monocle的一些关键步骤。

数据过滤（QC）

 1> HSMM<-monocle_cds
 2## 归一化 
 3> HSMM <- estimateSizeFactors(HSMM)
 4> HSMM <- estimateDispersions(HSMM)
 5#Filtering low-quality cells
 6> HSMM <- detectGenes(HSMM, min_expr = 3 )
 7> print(head(fData(HSMM)))
 8              gene_short_name num_cells_expressed
 9AL627309.1         AL627309.1                   9
10AP006222.2         AP006222.2                   3
11RP11-206L10.2   RP11-206L10.2                   5
12RP11-206L10.9   RP11-206L10.9                   3
13LINC00115           LINC00115                  18
14NOC2L                   NOC2L                 254
15> expressed_genes <- row.names(subset(fData(HSMM),
16+                                     num_cells_expressed >= 10))
17> print(head(pData(HSMM)))
18               orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters Size_Factor num_genes_expressed
19AAACATACAACCAC     pbmc3k       2419          779  3.0177759               1               1          NA                 779
20AAACATTGAGCTAC     pbmc3k       4903         1352  3.7935958               3               3          NA                1352
21AAACATTGATCAGC     pbmc3k       3147         1129  0.8897363               1               1          NA                1129
22AAACCGTGCTTCCG     pbmc3k       2639          960  1.7430845               2               2          NA                 960
23AAACCGTGTATGCG     pbmc3k        980          521  1.2244898               6               6          NA                 521
24AAACGCACTGGTAC     pbmc3k       2163          781  1.6643551               1               1          NA                 781

细胞分类（Classifying）

Monocle提供了一个简单的系统，可以根据您选择的marker基因的表达来标记细胞。您只需提供Monocle可以用来注释每个单元格的一组函数。

 1CDC20 <- row.names(subset(fData(HSMM), gene_short_name == "CDC20"))
 2GABPB2 <- row.names(subset(fData(HSMM),
 3                             gene_short_name == "GABPB2"))
 4
 5cth <- newCellTypeHierarchy()
 6cth <- addCellType(cth, "CDC20", classify_func =
 7                     function(x) { x[CDC20,] >= 1 })
 8cth <- addCellType(cth, "GABPB2", classify_func = function(x){ x[GABPB2,] >= 1 })
 9
10HSMM <- classifyCells(HSMM, cth, 0.1)
11table(pData(HSMM)$CellType)
12
13 CDC20  GABPB2 Unknown 
14      4      69    2565

我只是随意拿几个基因来试一下

1pie <- ggplot(pData(HSMM),
2              aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1)
3pie + coord_polar(theta = "y") +
4  theme(axis.title.x = element_blank(), axis.title.y = element_blank())
5

聚类

第一步是决定用哪些基因来进行聚类（特征选择）。我们可以使用所有的基因，但是我们将包括很多没有表达到足够高水平来提供有意义的信号的基因，它们只会给系统增加噪音。我们可以根据平均表达水平筛选基因，我们还可以选择细胞间异常变异的基因。这些基因往往对细胞状态有很高的信息含量。

1#Clustering cells without marker genes 
2
3disp_table <- dispersionTable(HSMM)
4unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
5HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
6plot_ordering_genes(HSMM)

红线表示单片基于这种关系对色散的期望。我们标记用于聚类的基因用黑点表示，其他的用灰点表示。

1# Plots the percentage of variance explained by the each component based on PCA from the normalized expression data using the same procedure used in reduceDimension function.
2# HSMM@auxClusteringData[["tSNE"]]$variance_explained <- NULL
3plot_pc_variance_explained(HSMM, return_all = F) # norm_method='log'
4

1> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 10,
2+                         reduction_method = 'tSNE', verbose = T)
3Remove noise by PCA ...
4Reduce dimension by tSNE ...
5> HSMM <- clusterCells(HSMM, num_clusters = 2)
6Distance cutoff calculated to 2.589424 
7> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
8+                    markers = c("CDC20", "GABPB2"))

我发现指定cluster为5的时候，只能画出来4个,为什么？

1HSMM <- clusterCells(HSMM, num_clusters = 5)
2plot_cell_clusters(HSMM, 1, 2)
3

Monocle允许我们减去“无趣的”变异源的影响，以减少它们对集群的影响。您可以使用到clusterCells和其他几个Monocle函数的residualmodelformula astr参数来实现这一点。此参数接受一个R模型公式字符串，该字符串指定要在cluster之前减去。

 1HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 2,
 2                        reduction_method = 'tSNE',
 3                        residualModelFormulaStr = "~Size_Factor + num_genes_expressed",
 4                        verbose = T)
 5HSMM <- clusterCells(HSMM, num_clusters = 5)
 6plot_cell_clusters(HSMM, 1, 2, color = "Cluster")  # 图不放了
 7
 8plot_cell_clusters(HSMM, 1, 2, color = "Cluster") +
 9  facet_wrap(~CellType)
10

monocle 支持半监督聚类Clustering cells using marker genes来定义细胞以及外部导入细胞类型（Imputing cell type）（cell与type对应文件），这里均不再展示。

构建轨迹

在发育过程中，细胞会对刺激做出反应，并在整个生命过程中，从一种功能“状态”转变为另一种功能“状态”。不同状态的细胞表达不同的基因，产生蛋白质和代谢物的动态重复序列，从而完成它们的工作。当细胞在不同状态间移动时，会经历一个转录重组的过程，其中一些基因被沉默，另一些基因被激活。这些瞬时状态通常很难描述，因为在更稳定的端点状态之间纯化细胞可能是困难的或不可能的。

单细胞RNA-Seq可以使您在不需要纯化的情况下看到这些状态。然而，要做到这一点，我们必须确定每个细胞的可能状态区间。

Monocle介绍了利用RNA-Seq进行单细胞轨迹分析的策略。Monocle不是通过实验将细胞纯化成离散状态，而是使用一种算法来学习每个细胞必须经历的基因表达变化序列，作为动态生物学过程的一部分。一旦它了解了基因表达变化的整体“轨迹”，Monocle就可以将每个细胞置于轨迹中的适当位置。

然后，您可以使用Monocle的微分分析工具包来查找在轨迹过程中受到调控的基因，如查找作为伪时间函数变化的基因。如果这个过程有多个结果，Monocle将重建一个“分支”轨迹。这些分支与细胞的“决策”相对应，Monocle提供了强大的工具来识别受它们影响的基因，并参与这些基因的形成以及如何进行分析分支。Monocle依靠一种叫做反向图嵌入的机器学习技术来构建单细胞轨迹。

在开始之前您也需要问，什么是拟时（伪时间）分析。

伪时间是衡量单个细胞在细胞分化等过程中取得了多大进展的指标。在许多生物学过程中，细胞并不是完全同步的。在细胞分化等过程的单细胞表达研究中，捕获的细胞在分化方面可能分布广泛。也就是说，在同一时间捕获的细胞群中，有些细胞可能已经很长时间了，而有些细胞甚至还没有开始这个过程。

当您想要了解在细胞从一种状态转换到另一种状态时所发生的调节更改的顺序时，这种异步性会产生主要问题。跟踪同时捕获的细胞间的表达可以产生对基因动力学一个大致的认识，该基因表达的明显变异性将非常高。Monocle根据每个cell在学习轨迹上的进展对其进行排序，从而缓解了由于异步而产生的问题。

Monocle不是跟踪表达式随时间变化的函数，而是跟踪沿轨迹变化的函数，我们称之为伪时间。伪时间是一个抽象的分化单位：它只是一个cell到轨迹起点的距离，沿着最短路径测量。轨迹的总长度是由细胞从起始状态移动到结束状态所经历的总转录变化量来定义的。

The ordering workflow

Step 1: choosing genes that define progress
Step 2: reducing the dimensionality of the data
Step 3: ordering the cells in pseudotime

基因的选择

首先，我们必须决定我们将使用哪些基因来定义细胞在肌生成过程中的进展。我们最终想要的是一组基因，能够随着我们研究过程的进展而增加(或减少)表达。

理想情况下，我们希望尽可能少地使用正在研究的系统生物学的先验知识。我们希望从数据中发现重要的排序基因，而不是依赖于文献和教科书，因为这可能会在排序中引入偏见。我们将从一种更简单的方法开始，但是我们通常推荐一种更复杂的方法，称为“dpFeature”。

1diff_test_res <- differentialGeneTest(HSMM[expressed_genes,],
2                                      fullModelFormulaStr = "~percent.mt")
3ordering_genes <- row.names (subset(diff_test_res, qval < 0.1)) ## 不要也写0.1 ，而是要写0.01。
4
5HSMM <- setOrderingFilter(HSMM, ordering_genes)
6plot_ordering_genes(HSMM)

根据时间点的差异分析来选择基因通常是非常有效的，但是如果我们没有时间序列数据呢?如果细胞在我们的生物过程中是异步分化的(通常是这样)，Monocle通常可以从同时捕获的单个种群重建它们的轨迹。下面是两种选择基因的方法，它们完全不需要知道实验的设计。

一旦我们有了一个用于排序的基因id列表，我们就需要在HSMM对象中设置它们，因为接下来的几个函数将依赖于它们。

降维

1#Trajectory step 2: reduce data dimensionality  
2HSMM <- reduceDimension(HSMM, max_components = 2,
3                            method = 'DDRTree')

排序

1#Trajectory step 3: order cells along the trajectory  
2HSMM <- orderCells(HSMM)
3plot_cell_trajectory(HSMM, color_by = "seurat_clusters")

有了树结构分类颜色是可以侄自己指定的。轨迹呈树形结构，Monocle事先不知道应该调用树的哪个轨迹中的哪个来调用“开始”，所以我们经常不得不再次使用root_state参数调用orderCells来指定开始。首先，我们绘制轨迹，这次按“状态”给细胞上色。

1plot_cell_trajectory(HSMM, color_by = "State")

1plot_cell_trajectory(HSMM, color_by = "Pseudotime")

“状态”只是单片树的术语，表示这段树。下面的函数便于识别包含时间为零的大多数细胞的状态。然后我们可以把它传递给orderCells

 1GM_state <- function(cds){
 2  if (length(unique(pData(cds)$State)) > 1){
 3    T0_counts <- table(pData(cds)$State, pData(cds)$seurat_clusters)[,"0"]
 4    return(as.numeric(names(T0_counts)[which
 5                                       (T0_counts == max(T0_counts))]))
 6  } else {
 7    return (1)
 8  }
 9}
10HSMM_1 <- orderCells(HSMM, root_state = GM_state(HSMM))
11plot_cell_trajectory(HSMM_1, color_by = "Pseudotime")

可以看出不同的排序条件拟时方向是不一样的。同时如果想看每个状态下细胞的分化可以指定分面，当然你也可以指定其他的分面方式。

1plot_cell_trajectory(HSMM_myo, color_by = "State") +
2    facet_wrap(~State, nrow = 1)

如果没有timeseries，可能需要根据某些标记基因的表达位置来设置根，使用您对系统的生物学知识。例如，在这个实验中，高度增殖的祖细胞群产生两种有丝分裂后的细胞。所以根应该有表达高水平增殖标记的细胞。我们可以使用抖动图找出对应于快速增殖的状态

1blast_genes <- row.names(subset(fData(HSMM),
2                                gene_short_name %in% c("CCNB2", "NOC2L", "CDC20")))
3plot_genes_jitter(HSMM[blast_genes,],
4                  grouping = "State",
5                  min_expr = 0.1)

单个基因的时间变化（我可以随意选择基因而你不可以，你要选有意义的生活）

1HSMM_expressed_genes <-  row.names(subset(fData(HSMM),
2                                          num_cells_expressed >= 10))
3HSMM_filtered <- HSMM[HSMM_expressed_genes,]
4my_genes <- row.names(subset(fData(HSMM_filtered),
5                             gene_short_name %in% c("YWHAB", "ICAM2", "ICAM2")))
6cds_subset <- HSMM_filtered[my_genes,]
7plot_genes_in_pseudotime(cds_subset, color_by = "seurat_clusters")